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Quasi-incompressible Multi-species Ionic Fluid Models 

Xiaogang Yangf Yuezheng Gong! Jun Li! Robert S. Eisenberg-and Qi Wang^ 


Abstract 

In traditional hydrodynamic theories for ionic fluids, conservation of the mass and linear momentum 
is not properly taken care of. In this paper, we develop hydrodynamic theories for a viscous, ionic fluid 
of N ionic species enforcing mass and momentum conservation as well as considering the size effect of 
the ionic particles. The theories developed are quasi-incompressible in that the mass-average velocity 
is no longer divergence-free whenever there exists variability in densities of the fluid components, and 
the models are dissipative. We present several ways to derive the transport equations for the ions, which 
lead to different rates of energy dissipation. The theories can be formulated in either number densities, 
volume fractions or mass densities of the ionic fluid components. We show that the theory with the 
Cahn-Hilliard transport equation for ionic species reduces to the classical Poisson-Nernst-Planck (PNP) 
model with the size effect for ionic fluids when the densities of the fluid components are equal and the 
entropy of the solvent is neglected. It further reduces to the PNP model when the size effect is neglected. 
A linear stability analysis of the model together with two of its limits, which is the extended PNP model 
(EPNP defined in the text) and the classical PNP model (CPNP) with the finite size effect, on a constant 
state and a comparison among the three models in ID space are presented to highlight the similarity and 
the departure of this model from the EPNP and the CPNP model. 
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1 Introduction 

Phase field models have been used successfully to study a variety of multiphasic phenomena like equi¬ 
librium shapes of vesicle membranes lfT3l [T4l . blends of polymeric liquids ll52l [53!. 54. 171 , multiphase 
fluid flows lfl9l l25l l34l l38l l35l l58l l57l l59l Ibll l63l. dentritic growth in solidification, microstructure evo¬ 
lution lf2Tl l40l [28l . grain growth J9J, crack propagation iflOl . morphological pattern formation in thin 
films and on surfaces [36, 451, self-assembly dynamics of two-phase monolayers on an elastic substrate 
ll37l . a wide variety of diffusive and diffusion-less solid-state phase transitions [H, 56], dislocation model¬ 
ing in microstructure, electro-migration and multiscale modeling ll49l . Multiple phase-field methods can 
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be devised to study multiphase materials |[57ll . Recently, phase field models are applied to study liq¬ 
uid crystal drop deformation in another fluid, liquid films, polymer nanocomposites, biofilms and cells 

HU m M [381 EO M EH [591 ED M (HI M M S3 MW1- 

Comparing to other mathematical and computational technologies available for studying multi-phase 
materials, the phase-field approach exhibits a clear advantage in its simplicity in model formulation, ease of 
numerical implementation, and the ability to explore essential interfacial physics at the interfacial regions 
etc. Computing the interface without explicitly tracking the interface is the most attractive numerical feature 
of this modeling and computational technology. Since the pioneering work of Cahn and Hilliard in the 50’s 
of the last century, the Cahn-Hilliard equation has been the foundation for various phase field models (7JIH1- 
It arises naturally as a model for multiphase material mixtures should the entropic and mixing energy of the 
mixture system be known. 

While modeling immiscible binary fluid mixtures using phase field theories, one commonly uses a la¬ 
beling or a phase variable (also known as a volume fraction or an order parameter) <|) to distinguish between 
distinct fluid phases. For instance c|> = 1 indicates one fluid phase while (f> = 0 denotes the other fluid phase in 
an immiscible binary mixture. The interfacial region is tracked by 0 < (|) < 1. For historical more than logical 
reasons, most mixing energies are calculated in terms of the volume fraction instead of the mass fraction in 
the literature lf20lfl2]l . Consequently, the system free energy including the entropic and mixing contribution 
has been formulated in terms of the volume fraction as well t20l fl21 . given in the form F[§, V(|), • • • ]. A 
transport equation for the volume fraction (]) along with the conservation equation of momentum and the 
continuity equation constitute the essential part of the governing system of hydrodynamic equations for the 
binary fluid mixture, where the volume fraction serves as an internal variable for the fluid mixture. 

In this formulation, the material incompressibility is often identified with the continuity equation 

V • v = 0. (1.1) 

This assumption is plausible and indeed consistent with the fluid incompressibility (11.11) only if the two fluid 
components in the mixture are either completely separated by phase boundaries when their densities are not 
equal or possibly mixed when the densities are identical. Otherwise, there is a potential inconsistency with 
the conservation of mass as well as conservation of linear momentum. This inconsistency has been identi¬ 
fied in lf38ll . but ignored by many practitioners using phase field modeling technologies for hydrodynamical 
systems. We note that this inconsistency occurs only in the mixing region of the two incompressible fluids, 
where the incompressibility condition ( 11.11) is no longer valid, indicating the mixture is no longer incom¬ 
pressible despite that each fluid component participating in mixing is. This type of fluids is referred to as 
quasi-incompressible in f38ll . A systematic fix to this problem for mixtures of incompressible viscous fluids 
was given by two of the authors in l3Tll . where the divergence free condition is modified to accommodate 
the quasi-incompressibility. 

In modeling of ionic fluids, one recognizes that the size of ions matters in most ionic solutions, in par¬ 
ticular in the ionic solutions in which life occurs, in the ocean, and of course in the very crowded conditions 
found in and near electrodes in batteries and electrochemical cells, in and around enzymes, ionic channels, 
transporters, and nucleic acids, both DNA and RNA lf68ll . Ionic solutions are hardly ever ideal: ionic size is 
almost always important. In multispecies ionic fluids above a certain concentration or under certain length 
scales, the size of the ions matters so that the same inconsistency issue in the models for ionic solutions 
arises again. That is one can not simply use the solenoidal condition in the velocity field as a proxy for the 
material incompressibility. A theory for multispecies ions of incompressible fluid flows that respects the 
material’s mass conservation and momentum conservation needs to be developed. 

This paper aims exactly at developing such a theory for a mixture of ionic fluid flows of multiple ionic 
species, in which the ionic densities are unmatched and different from that of the solvent, and their size 
effects are non-negligible. We require the theory to be dissipative while conserving mass and momentum. 
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One targeted application of this theory is in ion channel modeling Ifl5l[l6ll26ll . Ion channels provide enough 
data to distinguish between theories because measurements are available over a wide range of conditions 
E®. Hundreds of channel types are studied every day because of their biological and clinical significance 
1681 . Concentrations and electrical potentials are controlled in experiments and these provide sets of values 
for boundary conditions of mathematical models. Fitting the entire set with one set of structural parameters 
allows robust solutions of the inverse problem [5] 0 and thus allows models to be distinguished. Other 
applications of the model include electrolyte fluids, biological fluids with charged bio-species etc. This 
theory will be consistent with the mass and momentum conservation and demonstrates energy dissipation. 
In principle, a variety of transport equations can be developed for the ionic species should one knows the 
system’s energy dissipation rate. In this paper, we propose two types of transport equations based on a 
generalized Onsager principle lf60l . These two choices yield two types of species transport equations and 
corresponding energy dissipation rates. Their relations with respect to the existing electrolyte fluid models 
will be discussed in the text in details. 

The derivation follows the generalized Onsager principle approach OTl 60], leading to two types of 
transport equations for each ionic species in the form of Cahn-Hilliad and Allen-Cahn type equations, re¬ 
spectively. Apparently, these correspond to two distinct energy dissipation rates. Their applicability to real 
material systems can only be confirmed if one could measure the systems’ energy dissipation rates. How¬ 
ever, such measurements have not yet been made, as far as we know. So in most cases, people adopt one 
particular formulation over the others simply based on the leap of faith. 

For the new model, together with its limits in the extended Poisson-Nernst-Planck (EPNP) and the 
classical PNP with the size effect (CPNP), we will study their linearized stability on constant steady states. 
Instability of the PNP class of models is of direct biological interest. Actual biological channels invariably 
produce unstable currents iHdl that switch ’instantaneously’ between open and closed levels in a random 
telegraph process called single channel gating l24l . Instability in the models of this paper may turn into 
gating when the models are extended to include noise sources and are focused on the behavior of just one 
channel protein. However, we will not pursue the complicated issue in this paper; instead, we will focus 
on introducing the modeling framework and presenting a set of thermodynamically and hydrodynamically 
consistent theories, and discuss their predictions in a simple 1-D case to highlight the departure of several 
previously used PNP type models from the new model. 

The paper is organized as follows. First we present the mathematical formulation of hydrodynamic phase 
field theories for multispecies ionic fluid flows and various plausible formulations of the transport equations 
giving rise to the total energy dissipation. Then, we examine the theory in ID geometry to compare the 
theory with some existing PNP models with and without the size effect lfT5liT6t[26ll . Finally, we provide a 
concluding remark. 

2 Quasi-incompressible hydrodynamic models for ionic fluids 

We develop hydrodynamic models for a viscous, multispecies ionic fluid in an isothermal condition, in 
which mass, momentum conservation and the total free energy dissipation are preserved. The governing 
system of equations in the model includes the transport equations for all the ions, the Poisson equation for 
the electric potential, and the conservation equation for mass and linear momentum of the fluid, respectively. 

2.1 Mass and momentum conservation equations 

We first present the mass and momentum conservation equation. We consider the transport of viscous, 
ionic fluids made up of N different ionic species, each of which consists of a type of ionic particles of the 
identical size. Here, we tacitly assume the viscous solvent particle is a type of ions with a zero charge 
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ff30l[50 , 29J3). We denote the number density for each type of ions by /?,./' = 1, ■ ■ ■ ,N. The electric potential 
generated by these ionic particles is denoted by <J>. We denote the volume of each individual ionic particle 
by \’i and the mass by m t for i = 1, • • • ,N, respectively. Then, there is a constraint Y 4 L 1 n i v i = 1> which states 
that the excluded volume of the ions is a constant before and after the mixing. We identify i = a as the 
solvent component which is neutral. The total density of the mixture is defined by 

P = l£i mm. (2.1) 


We denote the intrinsic density of the zth species by p,- = which is a constant. Then, it follows that 


N 

P = L P««/Vi = Yj toPb 
1—1 


N 

L 

i—l 


( 2 . 2 ) 


where = tiiVi is the volume fraction of the ith ion in the mixture. We introduce the mass averaged velocity 
v. Then, the total mass and the linear momentum conservation yield 


| + V-(pv) = 0, 

P(^ + vVv) = V-t + fM 


(2.3) 


where T = — pol + l v is the total stress tensor, po is the hydrostatic pressure, i v is the extra stress tensor and 
is the interfacial force that yields the Ericksen stress for the mixture fluid system.We next turn to the 
derivation of the transport equations for the ions. 


2.2 Transport equations for the ions 

The free energy of the system is prescribed as F = J n f[n 1 , • • • . n,v]c/x, where fi is the material volume, 
and the density of the free energy functional is defined by [43 , 44] 


/[«!,••• ,n N ] =k B T'EHL 1 ^(lnzi,- l) + p e (4<J)„ + T> f ) + /^(x-y)G({?i;}^ 1 (x),{?i/}^ 1 (y))£/y, (2.4) 


where ks is the Boltzmann constant, T is the absolute temperature, N, is a generalized polymerization index 
for the zth ionic particle (N a = 1), p e = eo + Yj/Li Zien ,■ is the total charge density, z\ is the valence for type 
i ion and also denotes its sign (for solvent, we note that z a = 0), e is the unit charge, eo is the permanent 
charge density in the system, <!>„ is the electric potential generated by the total charge, <fy is a given external 
electric potential which is independent of the total charge and the total electric potential is <J> = d>„ + Ty. 
The first group in the sum represents the entropic contribution of the ionic particles to the free energy, the 
second part gives the electrical energy density of the system, and the third part gives the interaction of the 
excluded volume effect and the long-range interaction among the ions of finite sizes. 

The electrical energy density in the given external electric held is p e <fy and in the electric held generated 
by the charges is 2p r cp (| . The equations for the electric potentials Ty and <ty are 


f V • (eV<f>„) = — (<?o + ’LiZieni), ^ f V- (eV<D e ) = 0, 

1 < J > «lan = 0 1 1 ‘JMan = < £o(d^)> 


(2.5) 


where £ is the dielectric constant, <3>o is a given boundary function. Here the boundary condition is Dirichlet 
BC, it can be changed to other type boundary conditions. The external electric potential <fy is determined 
by the boundary condition with zero charge source. If <f>o = 0, there is no external electric potential. <fy is 
determined by the charge source with homogenous boundary condition and it can be expressed by using the 
Green’s function G(x,x') as 




( 2 . 6 ) 
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Then the variation of the electrical energy F e = f n p e (^<P n + & e )dx with the ion density n, is 

= z t e§ n + zie§ e = z,eT>. 

The equation for the total electric potential is 

f V • (eV<f>) = -(eo + 'LiZien i ), 

1 <fr|an = «fro@Q)- 

The third part of the free energy density can be approximated via expansions in a differential form 

/^(x-y)G({n,-}^ 1 (x),{n / }^ 1 (y))Jy«g[n 1 ;,;v.- ,n N ] =g({n,-}^ 1 ,{V/i,-}^ 1 ). (2.9) 

One specihc form of the function g accounting for the size effect of the ions is given by 

8 = k B T[ X5 =1 |n^ + Lti !l|Vn«|| 2 ], (2-10) 

where the coefficient matrix ^ is symmetry. The first part in the energy density represents a repulsive 
interaction due to the finite size effect while the second part is the conformation entropy associated with 
the heterogeneous distribution of the ions in space. This approximate function represents the lowest order 
approximation to the interaction potential with the long-range interaction, for which we will adopt in the 
rest of the paper. The chemical potential for the ith ionic particle is then given by 

Hi = H ( ,n n i) + £/^ ( /n/ - Y/'V 2 n,'] + ez,<S>. (2.11) 

Assuming there is no annihilation of charges between positive and negative ionic particles, each species’ 
charge and the total charge in the system is supposed to be conserved under the flux free boundary condition, 

f^n/dx = Cj, i = 1,••• ,N, Jq(LZi Zin t )dx = C = const, (2.12) 


where C,./ = I, • • ■ . /V and C are constants and C = 0 is called charge neutral. Indeed, annihilation can 
occur in biological systems and ordinary bulk ionic solutions when weak acids and bases (like acetic acid, 
i.e., vinegar, or sodium bicarbonate, i.e., baking soda) are involved as components of the solution or as 
side chains of the protein that forms the ion channel. Such effects arc significant in some cases, but they 
form a separate field of investigation, in theory, experiment, and indeed in medical practice, where they are 
particularly important. In this paper, we ignore those effects. 

We propose the transport equation for the ith ion as follows 

■ 

+ V • (vHf) = Bj,i = 1, - ,N , (2.13) 

where Bj is going to be determined from the total free energy dissipation in the following. We note that there 
arc two constraints of 5, as follows, due to the constraint of /;, and the total mass conservation, respectively. 
Using yOLi nm = Li^i ^ = 1> we have 


N N N 

i = 1 (=1 (=1 


It implies that 

V-v = lf =1 B/v, = Ef = i5,f. 
This gives us the first constraint on the B^s. 


(2.14) 


(2.15) 
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In addition, from the total mass conservation and p = m i n u we obtain 

N 

X>/tf/ = 0. (2.16) 

(=i 

This yields the second constraint on the 5-.v. The constraints warrants that the transport equations for each 
species are not completely independent. We next discuss two distinct ways to derive the transport equations 
for the ions and solvent following the generalized onsager principle ll60Tl . 

2.3 Formulation 1 

We denote the ath component (the solvent component) as the non-vanishing component in the mixture 
and then it follows from eq, (12.161) 


Ba — (2-17) 

The total free energy E = Jq(§ ||v|| 2 )£/x + F of the system consists of two parts: the kinetic energy and the 
Helmholtz free energy F. Now, we compute the total free energy dissipation rate as follows: 

§ = £/n[§ M 2 + f\<* 

= - /o[Vv : x - v ■ F« - £f =1 ^}dx + f dn n • (lf_, ^)dS 

(2.18) 

= - In [Vv : x - v • FW + i (V • vm + v • Vn ( ) - A nB,]dx 

= - /n( Vv : T v + l!Ll - £) -Hi + \Bi}dx, 

where d£l is the surface of the material volume Q, n is the unit external normal, the elastic force is identified 
as follows 


F( e ) 


N 

J> V *, 

i=l 


and the total pressure is given by 


N 

P = Po-Y.Ei’V- 

i=i 

In the last step, constraint eq. (12.171) is used. We also set the boundary condition 


n- 


5/ 

(Nn, 


= 0, 


(2.19) 


( 2 . 20 ) 


( 2 . 21 ) 


so that the surface integration is zero, i.e., Y!i=\ fdn n ' = 0- 

Next, we identify two forms of B, following the generalized Onsager principle to warrant energy dissi¬ 
pation of the system 11601 . They are associated with two famous transport equations: the Cahn-Hilliard and 
the Allen-Cahn equation, respectively. 
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2.3.1 Cahn-Hilliard dynamics 


In the first case, we choose Bj as follows 

= + for i*a, 


( 2 . 22 ) 


where the mobility coefficient matrix (kij. i. j / a) is symmetric and nonnegative definite. Then, using 
integration by parts, the energy dissipation rate is given by 


(2.23) 


f = - fa( Vy : T v + V[(-p) mi (± + %m a ] • hk 

V[{-p) m k(j- k -^)-W + ^Pa}}dx + surface term < 0 
provided Vv : x v > 0 and the surface term is zero. For viscous fluids, the viscous stress tensor is given by 


Tv = 2r| [D — \tr (D)I] + vr r (D)I. (2.24) 

where D = f (Vv +Vv r ) is the strain rate tensor, I is the identity tensor, q is the shear viscosity and v is the 
bulk viscosity. Then Vv : T v = 2r|D : D + (v — |r|)(fr(D)) 2 > 0 is satisfied so long as r| > 0 and v — |r| > 0. 
The zero surface term is warranted by the following no-flux boundary condition: 

n '{L*=i^V[(—£)-Aft + 5£ttx]} = 0. (2.25) 

We summarize the governing system of equations in this model in the following: 

|r + V.(vnO = -Li 1 V.^V[(-pK(^-X)_ w + ^ (x]5 for ,y a , 

V • v = -If i=1 Mj; ~ £)V • ^V[(-pK(i - ±)-p k + ^«], (2-26) 


P^F = V ■ [-(p ■+ I* t + Tv] + EjLi Ai/Vn,- = V • (-pi + x„) - Lf =1 zffV/if, 

and the equation for the electric potential is 

V • (eVffi) = -(e Q + Client). (2.27) 

where e is the dielectric constant. This model is not incompressible since V • v / 0 when densities are not 
identical. It is known as the quasi-incompressible model fTOll . This model is different from the previous 
models for ionic fluids. 

We remark that the previous models for ionic fluids assume the incompressible condition V • v = 0. This 
is valid only when p,- = p j,i,j= 1, •• ■ ,N. In this case, we end up with a self-consistent model as follows: 

^ + V• (vn,) = Lf = , v• X lk SJ\p k -Pal for i/a 


V • v = 0, 

pf = V -(-F I + Tv)-E£i« i -V rt , 

V- (eVffi) = -(eo + ZiZierii). 

In this model, the energy dissipation rate is given by 
clE f 

-TT = - / ( V v : T v + Y, V k’ - Pal ■ \kV[pk-P a]} < 0. 
dt iUa 


(2.28) 
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(2.29) 


For the above two model equation systems, the following boundary conditions are used: 
n 'W^i = 0 ’ 

(2.30) 

n • {It! h k n-p)Mj- k -£)-* + £**]} = 0. 

Together, they warrant that there is no boundary contribution to the energy dissipation and the constraints 
on the charge conservation in the system imposed by (12.12b are satisfied. The boundary condition for the 
electric potential is the Dirichlet boundary condition which is equal to a specified surface potential, and the 
boundary condition for the velocity held is the no slip boundary condition. 


2.3.2 Allen-Cahn dynamics 

Alternatively, we choose B, as follows 

B i = Zk=\hk\{-p)m k (j- k -±)-Vk + ! ^B a], for a, (2.31) 

where A* is the mobility coefficient, we obtain an Allen-Cahn type transport equation for the ith ion 

^ + V-(vn i ) = L^ =1 A*[(-p)OTt(^-^)-/ij t + ^/ia], for //a. (2.32) 

The other equations are given by 

V • v = -zu -£)-/*+ 


pf = V • l-(p + 1?=1 Wi) I + Tv] + If = i Ai/Vn/ = V • (-pi + Tv) - Lf =1 


(2.33) 


V-(eVffi) = -(eo + L iZietij). 

The boundary condition for this equation system is eq. (12.21b . The energy dissipation rate is given by the 
following 

f = -/n( v v:Tv+ 

1^=1 - l)-Bi + %ma]h k [(-p)m k (j- k -£)-* + %Ba]}dx < 0 , 

provided (A ; y) > 0. 

In the Allen-Cahn model, the charge conservation imposed by (12.12b may not be upheld. In order to 
impose the constraint approximately, we have to augment the free energy by adding a penalizing term 


N , 

!(/ ni-Q ) 2 
Jo. 


N 


Yzinidx — Cy 


(2.35) 


where L\ i are large positive numbers. An alternative approach is to enforce the constraints directly by using 
Lagrange multipliers in the free energy, 


Nr r N 

LiY( n i -C i )+L 2 { Y ZiUidx - c ), 
" Jo. \ 


(2.36) 


where L\ 2 are two Lagrange multipliers. These are common practices when one uses Allen-Cahn model to 
study multiphase fluid dynamics. We note that their physical validity is not widely accepted in the research 
community though. 

Note that Allen-Cahn and Cahn-Hilliard equations represent two different types of transport for scalar 
phase variables in a dissipative system |[39l . Higher order transport equations are also possible, but are rarely 
used. Thus, we will not pursue them in this study. 








2.4 Formulation 2 


By using constraint eq. (12.171) . we rewrite the energy dissipation rate as follows 
f = - /o{'W : Tv + If =1 K-p) % ~ MB<}dx 

= -/n{W : Tv + L^iK-^f -fii-Lmi]Bi}dx, 


(2.37) 


where L is a Lagrange multiplier, which is a function of the space and time. If we adopt the Cahn-Hilliard 
equation for the ionic species, the right hand term B, is chosen as 

Bi = ~ I^Li V • XijVK-p) ^ - q, - Lirij] , i = 1, - ,7V, (2.38) 

where A,,y is the mobility coefficient matrix. The constraint E/li m iBi = 0 implies 

N m 

£ V • ?1, 7 V[(-P)— ~ Bj — Lnij\nij = 0. (2.39) 

<J=i Pi 

It yields an elliptic equation for the Lagrange multiplier L: 

1/1,=! mmjV ■ XijVL = E-J =1 v. X tj V[(-p) g - //>/. (2.40) 


The Lagrange multiplier L is a solution of the elliptic equation. If the coefficient is a positive definite matrix, 
L is solvable in principle. In a special case where A,;/ are constants, the Poisson equation can be rewritten 
into 


V 2 L=[£ Xijnumj} 1 £ V-EyV[(-p)^ 

07= i i,j= i P/ 

Here, we don’t need to know the specific solution form for L. Then we have 

B, = ~ Lic=i V • hk{V[(-p)^~Bk] - ^ } 

= -Ef=! V-A*G*. 

The flux terms G^ are given by 

Gk = (-Vp)(— - ) - V/ft +-^= 1,2, ...,1V. 


(2.41) 


(2.42) 


(2.43) 


The terms wj = E) V _i L, 7 7r/,my, j = 1,2, ...,1V act as weighting factors. The difference between this model and 
the model derived in formulation 1 is that the correction factors are the weighted average terms. 

In a dilute solution, the solvent density is much larger than the other components, that is n a «/ for 
j / a. If we assume the mobility parameters A (/ - ~ where A,,- is a constant, then w ; - = Efli A/p/fyny ~ 

XjUjinj. Thus w a wy for j / a when mj and m a are not far apart, and this formulation reduces to the 
Cahn-Hilliard model derived in the previous subsection because 


l!J =1 Wjm k /pj 


~T7J - 

L"= i "7 


m k l!j= l Wjin^pj/mj 

Pa’ L“=i w/ 


'"a 


Vp 0 


G„ « 0. 


(2.44) 


For the solvent component, the governing equation of the density n a is 


+ V • (v/7 a ) — — V • ExaGoi ~ 0. 


(2.45) 
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Then we can drop the equation of the solvent component in our system and instead only consider the ionic 
components in this formulation. 

If we adopt the Allen-Cahn equation, the B, is chosen as follows 


B < = 1/ hj [{~P) 'jJ-Pj- Lm j \> 


(2.46) 


where A;y is the mobility coefficients. The constraint J^i-\ m i B i = 0 implies Lij=ihj[{-p)-p: ~ Pj ~ 
LmXm\ = 0. Thus, the Lagrange multiplier L can be solved as follows 


L = Ei,y=i hj mmj] Zij= 1 hi K-p) J- - Pj]m- 

The transport equation for the ith ion is given by 
a* V. V, r/ 


f+v-(vno =zti hk[(-p)(f- ^ 




Lj=i w i 


(2.47) 


(2.48) 


Using the same argument, if we assume the mobility parameters A y ~ A,/i/5,y, then vv ; = ^' v , A/yra/my ~ 
A jnjmj. Thus, w a » wy for j =4 a, which implies 


EjU ^ TOjt £^1 Wjin^j/mj ~ m 


Z.j=i w i 


Pa ’ 


Ty=i w t 


«<a 


Pa ■ 


and the governing equation of the solvent density /?„ is 


^ + V • (v/ia) « A aa [-p(^ - ^) -q a +q a ] = 0. 


(2.49) 


(2.50) 


dt 1 ’ • "u.u, l f^p„ p a 

This formulation reduces to the Allen-Cahn model derived in formulation 1. 

If (A ij) is a dense matrix, the two formulations are apparently different. However, if A,y = A8,-y, the 
Cahn-Hilliard equation derived in formulation 2 reduces to 

^T + V- (vn,-) = -AV 2 [—L/li - Pk + m i EfLi mmPi] • (2.51) 


p* ULi'A 

If m* = m,i= 1, • ■ - , N, it further reduces to 


^ + V • (m t ) = A V 2 [p k - 1 p,]. (2.52) 

Likewise, the Allen-Cahn equation reduces to 

^ + V • (\m) = -X\p k - jf Lf=i Pi}- (2.53) 

Both of these have been used by some researchers in the past to describe multiphase materials |f38l . 

Apparently, formulation 2 is different from formulation 1 and it seems to be a more general way of 
deriving the transport equations for the ionic species. However, if we choose L such that 


-p— -p a ~ Lm a = 0 (2.54) 

Pa 

and redefine 

1 N 

Ba = - V Bmu (2.55) 

m a trL 


we recover the model derived using formulation 1. This means that the transport equation for n a defined in 
reformulation 2 must be modified in order to recover the transport equation in formulation 1. However, this 
modification has no impact whatsoever on the energy dissipation rate. 

Another remark that we would like to make on these models is that each model yields an energy dissi¬ 
pation of its own. The choice of the model should therefore be made based on which energy dissipation rate 
best fits the real system to be modeled. 
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2.5 Model reformulation and reduction to existing models for multispecies ionic fluids 

The above models are formulated using number densities of the components in the fluid mixture. We 
can reformulate the model using the volume fraction <|),- or the mass fraction q since they are functions of 
the number density functions, (|), = n/V;, q = i = 1, • • • ,N, where v,- and //?, are constants, denoting the 
volume and the mass of each individual ionic particle, respectively. 

If pi = po,i = 1, • • • ,N, V • v = 0 and, in addition, we remove the entropic contribution of the solvent 
to the fluid mixture, i.e., we drop n a (\n n a — n a ), where a corresponds to the solvent component, from the 
free energy, the model reduces to the existing PNP model with the finite size effect l26l[27ll33ll . So, all the 
previous ionic fluid models can be regarded as the model applied to the case where all ions are of the same 
mass density and the solvent effect to the free energy is neglected. 

Next, we compare the new model with some of its limits and some existing models. 


3 Binary ionic fluid model 

We consider a mixture of two distinctive ionic components ( N = 3), where a = 3 corresponds to the 
solvent component, known as the binary ionic fluid model. The other two components in the fluid mixture 
are cations (positive ions) and anions (negative ions). We adopt the Cahn-Hilliard dynamics for the transport 
of ions. The governing system of equations is given by 

% + V • (v/ 1 /) = -v • hmn-pHij: - £)-Pi + fP3], 1 = 1,2, 

V. V = -- £)V• ^V[(-p)m,-(£ - £)-/* + §/i 3 ], 

(3.1) 

Pf = v-(-pi+^)-ELi« i v rt , 

V • (eVO) = -(e 0 + ZiZient), 

Here, we assume the mobility matrix is Xp = X,n,-8 (i , the mobility of each ion is only dependent on its own 
number density. The spatial gradients of the chemical potentials are given by 

Vqi = +eziVfl> + ksr[^i 2 V«2 + ^iiV«i — YiVV 2 «i], 

Vp 2 = ^^n2 + ez2^^ + kBT]^\2^n\ +^ 22^2 — 72 V V 2 n 2 ], (3.2) 


V » 3 = k B T ~ vlWn i~ V 2 Vn 2 . 

D 1 — V2«2 — Vi«i 

Where we assume that ^ 3 / = = 73 = 0, i.e., the interaction between the ions is dominant. The entropic 

contribution only shows up in the chemical potential of solvent (/q). 


3.1 Nondimensionalization 


We use a characteristic time scale to, length scale /o, and mass density scale po = P 3 , and the charac¬ 
teristic number density no to non-dimensionalize the physical variables. The mass density scale is chosen 
as the mass density of water here. Then, we denote the corresponding volume scale as vo = /q, mass scale 
mo = p 3 v’o■ The dimensionless variables are defined as follows: 


m 


nt_ 

no' 



h' 


to 

y= k, y ' p= 


p()/( 


2 Pi 


Pi = 


moli 


2 Pi- 


(3.3) 
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Then, the dimensionless parameters are given by 


1 V 3 5 1 

~ Tfl 

k B = -\k, 
mo /5 


H 13 1 


Po ’ 


Pl 

Po ’ 

li = 

-3* 

II 

U 5 * 

S3 

v = -^V 
P 0/5 

‘o e 

■*> = ^ 

I'll 

II 

p — m? 

moll 

t- - 9 9 

e%n 0 


V, 


(3.4) 


We set k B = 1 to obtain to = y -j—f and also set novo = 1 to obtain no = It’s easy to find that r'" = p,r]’ for 
i = 1,2 and rff = r v Q . For simplicity, we drop the ~ on the dimensionless variables and the parameters. The 
system of governing equations for the binary ionic fluid model in these dimensionless variables are given by 


^ + V • (vn;) = -V • hnFl-RiP ~ Vi + r>], * = 1,2, 

V • v = - lr =1 /?, V • hn.Vi-Rip - + rffi 3 \ = £ 2 =1 + V • (vn,)], 

pf = v -(-//! +Xv)-lL,n«Vp ( , 


(3.5) 


V • (eV<F) = ~(e 0 + £? =1 z,n,), 

where the parameters /?, = (^k)(4- — 1) for i = 1,2, the total mass density p = 1 — R\n\ —Rini, the solvent’s 
number density ns = r(j — r\n\ — r 2 n 2 . The spatial gradients of the chemical potentials are 

v A / i = jy^-Vni +ziV<J> + ^i 2 Vn2 + ^nVni -YiVV 2 m, 


= ^Vn 2 +Z 2 V<I> + ^i 2 Vni +^ 2 Vn 2 -Y 2 VV 2 n 2 , 


(3.6) 


V = -d' V »l-'W>n 

713 

In the following, we refer to the model as the full model, where the word ’’full” means that the model 
respects all conservation laws and accounts for the finite size effect and the solvent entropy. 


3.2 Models at regimes of two distinct length scales 

We examine the dimensionless full model at two distinct length scales. If we choose the length scale 
/(, = 10 9 m = 1 1 mi, we have the time scale t 0 = 1.55 x 10 M s. If we choose the length scale / 0 = 10 1 m = 
lOOnm, we have the time scale to = 1.55 x 10 ( \s. 

We set the first type ion is the positive ion with valence zi = +1 and polymerization index N] = 1; the 
second type ion is the negative ion with valence z 2 = — 1 and polymerization index /V 2 = 1. The values of 
the ratios of the ions’ volume, mass and density are tabulated in Table 1. The density ratio of the solvent 
and two ions is P 3 : pi : p 2 = 1 : 0.5 : 2, the volume ratio is V 3 : vi : v 2 = 1 : 2 : 1. The size differences of the 
three components are distinct. The parameters R\,R 2 are 0(10 2 ) in the smaller length scale /q = I nm. The 
compressibility of the flow (V • v / 0) in the full model can not be neglected. 

If the densities of the three components are the same, i.e. P 3 : pi : p 2 = 1 : 1 : 1, we have R\ = R 2 = 0, 
the flow becomes incompressible. Furthermore, when the density differences are distinct, but the larger 
characteristic length scale lo = 100 nm is used in the dimensionless system, the values of parameters R \. AS 
are very small, as in Table 1. If the corresponding terms of R , are dropped from the full model, the model 
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Table 1: The ratios of volume, mass and density 


Ratios 

Pi 

P2 

r\ 

r v 

r 2 

if 

rf 

Values 

0.5 

2 

2 

1 

1 

2 

Ratios 

r r> = ^0 

R 1 

Ri 

Values (/o = 1 nm) 
Values (Zq = lOOnm) 

40 

4x 10 7 

0.025 

2.5 x 10~ 8 

-0.025 
-2.5 x 10' 8 


reduces to a model that we call the extended PNP model (EPNP), in which the flow is incompressible: 

% + V • (v« ; ) = -v • XimVl-m + r>], i = 1,2, 

V • v = 0, 

(3.7) 

Pf = v -(-p I + Xv)-lLin i -V rt , 

V • (sV<f>) = -{e 0 + Y] =l Zitii). 

Here the mass density, the stress tensor and the chemical potentials are the same to those in the Full model. 

Furthermore, if we neglect the entropic contribution of the solvent to the fluid mixture, i.e., we drop 
n a (lnn a — 1), a = 3, from the free energy, we get the classical PNP model with the finite size effect (CPNP) 
lfl5l 1T6 L l22l [23l [5l l. This is equivalent to removing the /V 3 terms from the equations of the above EPNP 
model. To be clear, we note that the commonly used classical PNP model does not include the finite size 
effect. 

When the characteristic length scale used is /q = 1 nm, the parameters R 1 . Ri are not small. So, the full 
model must be used. The model is indeed different from the limiting PNP models even with the finite size 
effect. Note that this is the length scale regime that is applicable to the ion channel problem. The other 
model parameters are list in Table 2. 


Table 2: Model Parameters 


Symbol 

Parameter 

Value (Unit) 

Zo = 1 nm 

Zo = 100 nm 

0 

Shear viscosity 

1 x 10^ J kgm~ l s~ l 

15.54 

155.4 

V 

Bulk viscosity 

2.75 x 10~ 3 kgm~ l s~ l 

42.74 

427.4 

e 

Dielectric constant 

7.08 x 10 - w Fm~ l 

0.1145 

11.45 


Electric potential 

W 

38.65 

38.65 

Yi 

High order diffusion of 1th ion 

1.6606 x 10^“ / nv’mol- 1 

HP 1 

liP 

72 

High order diffusion of 2th ion 

1.6606 x 10- 27 m 5 mol~ x 

10~ 4 

10~ 14 

^1 

Mobility of 1th ion 

3.1083 x 10 11 kg~ l s 

0.02 

0.2 

^2 

Mobility of 2th ion 

3.1083 x 10 11 kg~ l s 

0.02 

0.2 


Self-interaction of 1th ion 

1.6606 x 10~ 5 

1 

10- 6 

^22 

Self-interaction of 2th ion 

1.6606 x 10- 5 rrv'mol _1 

1 

10^ 6 

in 

Interaction of the two ions 

1.6606 x 10- 5 nr’mol -1 

1 

10~ 6 


3.3 Comparison of the full model with the limiting PNP models in ID space 

We compare the full model with the EPNP and CPNP models in ID space, assuming the system is 
homogeneous in the (y,z) directions and depends only on x and time t (i.e., the variables are functions of 
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(, t,x ).) The domain for x is assumed finite given by = [0,L r ]. The governing equations of full model in 
ID are given explicitly by 

^ + MO' = -{Mi Mi (/>)' - (w)' + 'T(w)']}', 

^+(yi/12)' = -{^n 2 [-R2(py - oey+^cusm', 

(vi)' = ^[^ + (vim)']+^[^ + (vi« 2 )'], 

p(^ + v,(vl)') = -(p)' + (fl^+v)(v 1 )"-[ ?^ l( ; Ul)' + « 2 ( ; u 2 )' + n3(A'3) , ], (3-8) 

p(^ + v,(v 2 )')=tl(v 2 )", 

P(^ + Vi(v 3 )')=tl(v 3 )", 

(<P)" = -[e 0 + E? =1 z;«i]/e, 

where (•)' = (•)" = , and the gradients of the chemical potentials are 

(/rO^^^iy+zi^y+^^y+^iiCniy-yi^i)''', 

{&)'= ^(«2)' + Z2( < S > )' + ^12(«l)' + ^22(«2)'-Y2(«2) W , (3.9) 

( w y = lMM. 

The unknowns are ni,« 2 ,p,vi,v 2 ,v 3 ,‘I>, which are fully coupled. 

The ID EPNP model is much simpler now. First, from the incompressible condition (vi)' = 0, we hnd 
that vi = 0 due to the fixed boundary condition of v. Then, the pressure p and v 2 ,v 3 are determined from 
the momentum equation. The independent unknowns in the EPNP model are then «i,/r 2 ,<t> and the ID 
governing equations are given by 

^=-{*.i»i[-(m)'+>r (»)']}', 

^ = -{Wi2[-Cu2)' + <(»)']}', (3-10) 

where the gradients of the chemical potentials are given by (13.91) . 

If we further remove the p 3 terms from the ID governing equations of the EPNP model, we get the ID 
equations of the classical PNP model with the finite size effect (CPNP). In the following, we will compare 
these three models explicitly in ID. First, we examine their linear stability properties. 

3.4 Linear stability of the constant state 

If we assume eo = 0 (namely, the system does not have a permanent charge present), there exists a 
constant solution of the full model, which is a solution of all limiting models, n\ = hi = n°, v = 0 , p = 
0, = 0, where n° is constant such that 

= r ' 0 - r\n° - r' 2 n° > 0. (3.11) 
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This inequality is necessary to ensure that the solvent density is greater than zero. We perturb this constant 
solution as follows: 


n i = n° + £e at+ikx n ( l , n 2 = n° + Ze M+ikx n ° 2 , 

V| = gg«f+ifacy0 p _ £ e ca+ikXpO q, _ ggCW+i'fctqjO 


(3.12) 


Here e <C 1 is a small parameter, a is the growth rate and k is the wave number. First, we point out that the 
velocity components \ 2 , V 3 are decoupled from the rest of the system in the linearized equations and they do 
not contribute instability in this problem; so, we only consider the coupled system involving the remaining 
variables: The linearized eigenvalue problem for the Full model is given in the Appendix. 

The asymptotical analysis in the small wave number regime shows that the instability can incur only when 
£12 is negative enough. But, > 0 in the model. So this mode of instability is absent from the full model 
and its limits. The system is stable for long wave (small wave number) perturbation. It is easy to find that 
the system is also stable for short wave (large wave number) perturbation. From the numerical studies, we 
find that the intermediate wave instability appeal's when £,12 is sufficiently large. The analytical result of the 
intermediate wave instability is hard to obtain from the full model, but easy from the EPNP model. We thus 
focus on the linear stability of the limiting EPNP model in the following. 

The linearized eigenvalue problem is given in the Appendix. The instability condition is A < 0, where 


A — A-i X 2 (n 0 ) 2 1 1 ak 2 + bk 4 + + £11 + -^r) + Yt( 


1 


+ ^22 + ^)]A: 6 +YiY2^}<0. (3.13) 


Here, the parameters a, b are defined by 

= [jvJp + + ^22 + 2^12 + 

" ™ ,, r Vy.m 

+ fe2 + Tr)w7T+ 


a ■ 


u — Y1+Y2 , 1 , re,, , ’Vi \ _J_ 

e “T" NiN 2 (n u ) 2 > N 2 n' 

ry^22+r;^n , e e 
-^-+ 41^22 -S 12 -—^ 


,yi‘ 

d 

2 " 


s!2- 


Because A,i, A, 2 ,Yt»Y 2 ,^ 11,^22 are a ll positive, so the coefficients of k s ,k 6 are all positive. It implies that 
A > 0 for large wave numbers. This means that the system is stable for short waves. For long waves (small 
wave numbers), since the parameter a > 0, then A > 0 for small \k\. So, the solution is stable. Analogously, 
the mode of instability is absent from the CPNP model, which is a limit of the EPNP model, at both long 
and short waves. 

For intermediate waves, we notice a possible instability if b is negative, i.e., 


a > 0 , b < 0 . 


(3.14) 


In certain parameter regimes, the growth rate oq (given in the Appendix) can be negative for a very small \k\, 
becomes positive for some intermediate values of \k\, and then turns to negative again at large \k\. Assuming 
Yi = Y 2 = § 1, we obtain the roots of A = 0 asymptotical. Then, we obtain the cutoff wave numbers. We 


denotec =^ + ^11 + ^ + 


A^U+^22 + : 


, then we have 


A/(kiX 2 (n°) 2 k 2 ) = a + bk 2 + cbk 4 + $ 2 k 6 . 


There are two positive roots of k z , corresponding to two cutoff wave numbers k\% o11 , asymptotically: 

K 0ff ) 2 = -f -#8 + 0(5 2 ), (k c 2 utoff ) 2 = x i+ Xl +0( 5), (3.15) 


where xq = 


_ -c+yc 


- 4 b 


,X\ = 


b+ 2 cxo+ 3 xk 


. We only retain the first two terms in the asymptotic roots. The 


parameter A is negative when the wave number is between the two cutoff wave numbers 0 < k x 


cutoff 


<k< 
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k™ to ff' The growth rate ai is positive in this intermediate wave number regime. This instability depends 
strongly on the interaction parameter ^ 12 , the instability condition is satisfied for a sufficiently large £, 12 . 
We fix the parameter N\ =N 2 = l,A,i = X 2 = 0.02,yi =''{ 2 = 10 4 .4 11 = ^22 = 1,/?° = 1 in the length scale 
/o = I/7/77 regime, and vary the parameter £ 12 . When ^12 > 2.06, b < 0, the intermediate wave instability 
incurs. In FigureQJa), we show the curves of the two cutoff wave numbers as a function of ^ 12 . The smaller 
cutoff wave number k\ uto ^ is decreasing and the larger cutoff wave number k(“" >! ' is increasing as £12 
increases from 2.046. The unstable wave number regime (^ uf °//^™/o//) w j ( j ens as q |2 increases. 




Figure 1: Cutoff wave numbers as functions of £,12 with n° = 1 and as functions of n° with £,12 = 2.2. 


When 6 —>• 0, we have k\”'"^ —> y/—f and k) 1 " 0 ^ 1 —> + 00 . The unstable wave regime is k > k c “ to ^. 
The system is unstable for large wave numbers (short waves), which is known as the Hadamard instability. 
Hence, the high order diffusion coefficients yi, y? have the effect to suppress the short wave instability. 

This intermediate wave instability is also dependent of the constant state n°. When the interaction 
parameter ^12 = 2.2 is fixed, but n° is varying, we find that b is positive for small n° and negative for large 
n°. That means the system is stable for dilute solution but unstable for rich solution. We also plot the cutoff 
wave numbers as functions of n° with fixed ^12 = 2.2 in FigureQJb). The instability appears when n° > 0.87, 
and the unstable wave number regime widens as n° increases. 

This intermediate wave instability is a feature of these three models. Through a numerical investigation, 
we confirm that this instability property can occur in all three models. In the following example (Figure 
121), we use parameter values N\ = N 2 = l,A,i = = 0.02,yi = y> = 10~ 4 ,^n = ^22 = 1,^12 = 2.2,/?° = 1 

in the length scale Iq = 1 nm regime. The instability condition (13.141) is satisfied. The two asymptotical 
cutoff wave numbers of EPNP model are k c “ to ^ = 9.24 and / = 45.23, respectively, when length scale 
/o = In///. For the three models, the relation between the length scale and the growth rate follows a simple 
scaling law: we denote the two length scales as the corresponding growth rates as and 


A2) 

the cutoff wave numbers as respectively. If = K, then the cutoff wave number ratio follows 


F 2 > _ 


K while the growth rate ratio follows 


a f\Kk) 

a i 11 (*) 


= K 2 - 5 . 


This can be inferred from the definition of time 


scale t 0 = ^ 


,2.5 

l o 


, (mo Iq). The numerical results in Figure [2] also confirm this analysis. 

The analysis and numerical results show that the growth rates can be positive in some intermediate wave 
number regime depicted in Figure [2l instead of near the zero wave number range. In this case, the growth 
rate of the Full model is the smallest while the EPNP model’s is the highest. From the linear stability 
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analysis, we notice that this instability is associated with a large interaction parameter ^ 12 , a consequence 
of the finite size effect. A positive ^12 means that the interaction between different species due to their 
steric effects is repulsive. The analysis and numerical results tell us that the intermediate wave instability 
appears when the repulsive effect is sufficiently strong in the three models. This also can be obtained from 
the interaction free energy density g. The repulsive interaction due to the finite size effect is represented by 

+2^i2«i«2 + ^22?It)’ which can be rewritten as + -^=m) 2 + (£,22 — )n\ )■ When 

£12 is sufficiently large, c, 11 £,22 < 0, this quadratic form is hyperbolic type without lower bound. In 

the next nonlinear simulations, we only consider the cases ^ 11^22 > 0 , with out the intermediate wave 

instability. 





Figure 2: The growth rates of the full model and the two PNP models with length scale /q = 1,10, 100«m, 
respectively, in the parameter regime of intermediate wave number instability. The values of growth rates 
increase as the length scale /(> increases. The two cutoff wave numbers of the EPNP model in the length 
scale Iq = him regime are 9.24 and 45.23, respectively. The full model is more stable than the other two 
models in this regime. 


3.4.1 Discussion on the finite size effect 


The hard sphere repulsion characterizes the finite-size effect of ions, witch keeps ions apart. The free 
energy density due to the finite-size effect is 


fK(x - y)G(Mf =1 (x),{ir,}f =1 (y))dy = f £f =1 


2 |x—y| 12 


-nAxin 


(y)dy, 


(3.16) 


where n, and aj are the radii of ion i and j, and e,-j is the energy coupling constant between ion i and j. Thus, 
in the free energy function, we have the convolution integral with the following form 

11 ^z^ni(^)nj(y)dydx. (3.17) 


We can approximate the above convolution integral by truncating the kernel , with the cutoff length 

1^ y I 

8 . As discussed in the paper (69J, when the cutoff length 8 goes to zero, this convolution integral can be 
approximated by the integral 


S & fni(x)nj(y)dx, (3.18) 

with .S '5 ~ 8 121,1 1 where d is the dimension. The free energy density due to the finite-size effect can be 
written as 

E/j =1 ^(«i + aj) 12 S$ni(x)nj(x) = k B T^ J=l \mnj, (3.19) 
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with q,y = jjjEjjicii + a ; -) 12 5§. We add the conformational entropy in terms of the derivative form to com¬ 
pensate for the approximation error, then the energy density for the finite-size effect is approximated by 

8 = WLf/ =1 fywj + Hi !l|Vn,-|| 2 ], (3.20) 

where y,■ is a small parameter, witch can be zero. In the paper ll69l , the following £,y values for the cross 
hard-sphere potential terms for some familiar ions (Na + ,C1~,Ca 2+ ) are used: 

£-Na,Na ■ £ cl,cl '■ £-Ca,Ca ■ £-Na,ci '■ £-Na,Ca ■ £ci,Ca = 1:1:1: 0.955 : 1 : 0.961. (3.21) 

Also in the paper lf69l . the ratios of the interaction coefficients ^yy are given for some familiar ions (Na + .Cl .Ca 2 ] ) 
as follows 


^Na,Na '■ 4 Cl,Cl '■ 4 Ca,Ca '■ 4 Na,Cl ■ 4 ,Na,Ca ■ 4 Cl,Ca = 1 : 2280 : 1.64 : 42.2 : 0.642 : 50.4. (3.22) 

It is easy to verify that £, 11^22 > 0 for two of the three ions. For the familiar ions, the interaction 

coefficients y are in the stable regime. That is the reason we only consider the stable cases in the nonlinear 
simulations next. 

3.5 Nonlinear dynamics 

We next explore nonlinear dynamics of the models in the linearly stable regime. We use the characteristic 
length scale /o = I run and set the domain as x E [0,10]. The values of the interaction parameters are chosen 
as = ^22 = l,4i2 = 0-8, satisfying ^ 11^22 — 4p > 0. We also set diffusion coefficients Yi = 72 = 0. 
The boundary conditions for the number densities n\ ,ti 2 are no-flux boundary conditions (12.251) : for the 
velocity V|, the boundary conditions are set at vi | x =o.io = 0, and for the electric potential <E>, they are set at 
ffi| x= o = 0, ffi 1.^,,) = ffi 0 , where <J>o is the electric potential at the right boundary x = 10. We set <f>o = 1 in the 
following simulations. The initial conditions are given by n\ = m = n° = 1, Vi = 0, p = 0, <I> = ffio.r/10. 
The given external electric potential is dy. = ffio.r/10. The dimensionless mobilities are given as \\ = Xi = 
0.02. We compute the ionic number densities using the Full model, the EPNP model and the CPNP model, 
respectively. 

Figures [3] depicts the final steady states of the Full model and the EPNP model, and the difference 
between them, where = ^22 = 1-4)2 = 0.8 in the stable regime. The states of the number densities are 
almost identical in the middle of the domain, while the visible differences appear near the two boundaries. 
Because the electric potential is positive at the right boundary and zero at the left boundary, some negative 
ions gather at the right side while positive ions gather at the left side due to the Coulomb force, forming two 
visible boundary layers. 

As shown in Figure [3] the density differences between the two models are about 0(1) x 10 2 near the 
boundaries. As a conclusion, the compressibility of the flow in the full model plays relatively important 
role, it impacts the aggregation effect of the ions near the boundaries. 

In the above example, the density ratio is chosen as P 3 : pi : p 2 = 1 : 0.5 : 2. By halving density pi and 
doubling density P 2 to increase the density differences, we reset the density ratio as P 3 : pi : P 2 = 1 : 0.25 : 4 
and P 3 : pi : P 2 = 1 : 0.125 : 8, while maintaining the volume ratio unchanged at V 3 : vi : V 2 = 1 : 2 : 1, then 
the dimensionless parameters are A'| = 0.0375. AS = —0.075 and A’| = 0.04375. AS = —0.175, respectively. 
In these cases, the size differences of the three components become larger. As shown in Figure |4j the 
differences between the Full model and the EPNP model become larger as the size differences become 
larger. When we halve the density pi and double the density P 2 , the absolute maximum difference between 
the Full model and the EPNP model is almost doubled. As a result, when the size differences between the 
components are enlarged, the parameters R] M 2 are no longer small so that the compressibility of the flow 
can no longer be neglected. 
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(a) Number density (b) Differences between the 
models 




(c) Electric potential <t> 


(d) Energy F(t) — F( 0) 


Figure 3: Steady states of the ionic densities and the electric potential of the Full and the EPNP models with 
= £,22 = 1,^12 = 0.8 in the stable regime, respectively. The differences appear near the boundary with 
absolute maximum difference 0.011. The curve of energy difference F(t) —F( 0) is plotted with respect to 
time. The total free energy F(t) decays to a constant when the final steady state is obtained. 



XXX 


(a) p 3 : pi : p 2 = 1 : 0.5 : 2 (b) p 3 : pi : p 2 = 1 : 0.25 : 4 (c) p 3 : p, : p 2 = 1 : 0.125 : 8 

Figure 4: Differences between the Full and the EPNP model with = ^ 22 = 1,^12 = 0.8 in the stable 
regime, (a) The density ratio is p 3 : pi : p 2 = 1 : 0.5 : 2 and R\ = 0.025. AS = —0.025, the absolute maximum 
difference is about 0.011; (b) the density ratio is p 3 : pi : p 2 = 1 : 0.25 : 4 and R\ = 0.0375,f? 2 = —0.075, 
the absolute maximum difference is about 0.022; (c) the density ratio is p 3 : pi : p 2 = 1 : 0.125 : 8 and 
R\ = 0.04375. AS = —0.175, the absolute maximum difference is about 0.044. 


We also compare the EPNP and the CPNP model in the stable regime with = ^ 22 = 1,^12 = 0-8 

in Figure [5] The density ratio p 3 : pi : p 2 = 1 : 0.5 : 2 is used. The values of the parameters are set at 
r'" = I, > J 2 = 2. The number density differences between the two models are about 0(1) x 10 2 near the 
boundaries. The differences near the boundaries in n 2 are bigger than that in n\ . The reason is that in 
the CPNP model, the term r^/u 3 is dropped in the n, transport equation. In this example, P,” > rf. so 
the differences near the boundaries of n 2 are bigger. Consequently, the solvent’s chemical potential q 3 in 
the EPNP model plays a relatively important role, it impacts the aggregation effect of the ions near the 
boundaries. 

Next, we consider ionic concentrations without the finite size effect and compare them with ionic con¬ 
centrations with the finite size effect using the classical PNP model. In the following, the OPNP means the 
classical PNP model without finite size effects (i.e., £,n = q 22 = £12 = 0 and Yi = y 2 = 0.) The differences 
also appear in the areas near the two boundaries. As shown in Figures [6] the differences of the ionic density 
can reach up to 0(1) x 10 1 . The finite size effect plays an important role in the system, it impacts the 
aggregation effect of the ions near the boundaries, as studied in the papers If26ll27i [33l [691 . 

Based on our numerical investigations and the linear analysis, we conclude that the ID steady states 
of the number densities are nearly identical in the middle of the domain in all three models in the stable 
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X XX 


(a) Number density (b) Difference between the models (c) Electric potential <f> 

Figure 5: Steady states of the ionic densities and electric potential of the EPNP and the CPNP model with 
^11 = ^22 = I.^i 2 = 0.8 in the stable regime. The difference appears at the boundary layers near both ends 
of the domain. 



X XX 


(a) Number density (b) Difference between the models (c) Electric potential <t> 

Figure 6: Steady states of ionic densities and the electric potential of the CPNP and the OPNP model. In the 
CPNP model £,n = ^22 = I. c, 12 = 0.8 in the stable regime. The difference appears at the boundary layers 
near both ends of the domain. 

regime. The differences lie in the areas near the boundaries. The compressibility of the flow, the chemical 
potential of the solvent and the finite size effect are three main reasons that lead to the differences. So, our 
quasi-incompressible model (the full model) seems to be more reasonable because the mass and momentum 
conservation laws are preserved in the model while the other models don't respect the two fundamental 
physical conservation laws. 

Further investigations in higher dimensions is necessary to evaluate the difference among the models, 
which will be conducted in a sequel. 

4 Conclusion 

We have developed systematically a set of quasi-incompressible theories for ionic fluids of multiple 
species that respect not only momentum conservation but also mass conservation at the presence of the 
ionic species. The previous PNP type models are approximations of the more fundamental theories when 
densities of different ionic species are distinct. In these theories, we consider the entropic contribution from 
each ionic species together with the ion-ion interaction due to the finite size effect. The limiting cases include 
the extended PNP, the classical PNP with the finite size effect, and the classical PNP model without the finite 
size effect. At the length scale larger than hundreds of nanometers, all models agree with the classical PNP 
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model very well. At the length scale in a few nanometers, the models can predict quite different stability 
behavior for homogeneous equilibrium states. In nonlinear dynamics, the ionic number densities are nearly 
identical in the middle of the domain, but the differences lie in the areas near the boundaries. Apparently, 
three main factors in the compressibility of the flow, the chemical potential of the solvent and the finite size 
effect of the ions can lead to the discrepancy in model predictions. We tend to believe that the new model is 
more accurate since it obeys the two fundamental physical conservation laws in mass and linear momentum 
while the others don’t. 
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5 Appendix 


The linearized eigenvalue problem for the Full model is formulated as follows, 
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where the parameter values are given by n ( \ = — r\n° — r'^rP 

nents in the matrix are defined as follows 
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(5.1) 


— R 2 n . The other compo¬ 


rt = —ikr\ + ikn°(j±p + + £12) + /k 3 n°Yi, 

B 2 = —ikr\ + ikn°(j±p + 1, 22 + £12) + ik 3 n° y 2 , 

Cl = XjuV/Vj’ + Xi n° (^ + £11 )k 2 + Xi«°Yik 4 , 

C 2 = Xin°r^rl^ J + 'kin 0 (Z 3 n)k 2 , 

D\ = X2n 0 r^r v 1 ^j + X 2 n 0 (l ) i2)k 2 , 

Do = WVTr!;|j + X 2 n°(jj^p + £ 22 )k 2 + l 2 n°y 2 k 4 . 


Although the coefficient matrix is 5 x 5, the characteristic polynomial of the coefficient matrix is a third order 
polynomial of growth rate a, which yields three independent eigen-modes. Using an asymptotic analysis at 
small wave numbers \k\ <C 1, the three asymptotic growth rates are obtained asymptotically: 


ai = T\k 2 + Oik 4 ), 

aM = _ wp; + 0 ( t2 ), 


(53) 
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where 


Tl ~ “(Xj+CM [ p £"3(^7? + + ^n +^22 + 2^12) + Pa(?'! + r ' 2 )(( f J i + ^) + (^1 +^2)^3)], 

T 2 = A,i?l2Pi:(« 0 ) 2 (^l +i?2) 2 + £; 

T 3 = -4(n°) 2 p k £{X l +l 2 )(l l R 2 l +X 2 Rl), 

T 4 = 2n°p k £(l l R 2 + X 2 Rl). 

Notice that 0 * 2,3 < 0 for small k due to 7) < 0. The eigenvalue oq > 0 when T\ > 0, 
a ^ + a ^ + ^1+^22 + 2^12 <+ r^) + {Ri+R 2 )n Q 3 ). 


(5.4) 


P/S»3 


(5.5) 


This is the instability condition for long waves for the Full model. It follows from eqn (13.111) that (V" + 
t ~2 ) + (R 1 +R 2 )n 3 > 0. So, the instability can incur only when £,n + £,22 + 2c, 12 is negative enough. But, 
q,-y > 0 in the model. So this mode of instability is absent from the full model. 

For the EPNP model, only n\ ,n 2 ,<t> are coupled, the eigenvalue problem is given by 


0 0 0 \ / -k 2 i -i 

0 a 0 + n°k 2 Ci C 2 

0 0a) I -h.n°k 2 D ] D 2 
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= 0 , 


where 


Ci = Xin°r^r\^+hn°(j^ + ^n)k 2 + 'kin 0 Yik 4 , C 2 = +lin°t, l2 k 2 , 

D\ = X 2 n°r%r\^ +’k 2 n%\ 2 k 2 , D 2 = l 2 n Q ^r^ + l 2 n 0 (j^+^ 2 2)k 2 + l 2 >i°y 2 k 4 
Eliminating <f>°, the system reduces to 


a 0 
0 a 
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D\ - X 2 rf\ D 2 + X 2 n' 


oi 
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,0 1 


= 0 , 


The characteristic polynomial is quadratic and given by 

a~ + [Ci + D 2 + (Xi + A,2)n^f]oc+A = 0, 

A = CiD 2 - C 2 D x + Pl 2 (Ci +C 2 ) + Xi(Di +D 2 )]n°l 
The two growth rates are given by 


oq = —2A 


(Ci + D 2 + (Xi + A, 2 )n°f) + sj (Ci +D 2 + (A.i + X 2 )n°^) 2 — 4A 


a 2 = 


(Ci +D 2 + (Xi +A,2)»°p) + \ (C\+D 2 + (^-1 +^2)* 7 °p) 2 — 4A 


(5.6) 


(5.7) 


(5.8) 


(5.9) 


(5.10) 


Notice that Ci >0,D 2 > 0. So, Re( a 2 ) < 0 and Re{ ai) can be positive only if A < 0. 
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